Libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('keep.trailing.zeros',TRUE)
Reporting the Metrics
ROCAUC <- NULL
CstatCI <- NULL
LogRangp <- NULL
Sensitivity <- NULL
Specificity <- NULL
for (topf in topFive)
{
CstatCI <- rbind(CstatCI,RRanalysis[[topf]]$c.index$cstatCI)
LogRangp <- rbind(LogRangp,RRanalysis[[topf]]$surdif$pvalue)
Sensitivity <- rbind(Sensitivity,RRanalysis[[topf]]$ROCAnalysis$sensitivity)
Specificity <- rbind(Specificity,RRanalysis[[topf]]$ROCAnalysis$specificity)
ROCAUC <- rbind(ROCAUC,RRanalysis[[topf]]$ROCAnalysis$aucs)
}
rownames(CstatCI) <- topFive
rownames(LogRangp) <- topFive
rownames(Sensitivity) <- topFive
rownames(Specificity) <- topFive
rownames(ROCAUC) <- topFive
pander::pander(ROCAUC)
| age |
0.589 |
0.491 |
0.686 |
| wt.loss |
0.560 |
0.462 |
0.659 |
pander::pander(CstatCI)
| age |
0.559 |
0.560 |
0.499 |
0.619 |
| wt.loss |
0.515 |
0.515 |
0.456 |
0.574 |
pander::pander(LogRangp)
pander::pander(Sensitivity)
| age |
0.1157 |
0.0647 |
0.187 |
| wt.loss |
0.0496 |
0.0184 |
0.105 |
pander::pander(Specificity)
| age |
0.872 |
0.743 |
0.952 |
| wt.loss |
0.894 |
0.769 |
0.965 |
meanMatrix <- cbind(ROCAUC[,1],CstatCI[,1],Sensitivity[,1],Specificity[,1])
colnames(meanMatrix) <- c("ROCAUC","C-Stat","Sen","Spe")
pander::pander(meanMatrix)
| age |
0.589 |
0.559 |
0.1157 |
0.872 |
| wt.loss |
0.560 |
0.515 |
0.0496 |
0.894 |
Cox Model Performance
Here we evaluate the model using the RRPlot() function.
The evaluation of the raw Cox model with RRPlot()
Here we will use the predicted event probability assuming a baseline
hazard for events withing 5 years
timeinterval <- 2*mean(subset(lung,status==1)$time)
h0 <- sum(lung$status & lung$time <= timeinterval)
h0 <- h0/sum((lung$time > timeinterval) | (lung$status==1))
pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
Initial Parameters
| 0.85 |
578 |
index <- predict(ml,lung)
rdata <- cbind(lung$status,ppoisGzero(index,h0))
rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=lung$time,
title="Raw Train: Lung Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






By Risk Categories
obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
121 |
100.4 |
144.6 |
77.3 |
1.57 |
1.30 |
1.87 |
3.70e-06 |
| low |
83 |
66.1 |
102.9 |
51.6 |
1.61 |
1.28 |
1.99 |
4.97e-05 |
| 90% |
38 |
26.9 |
52.2 |
25.9 |
1.47 |
1.04 |
2.01 |
2.33e-02 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event Analysis
isevent <- rrAnalysisTrain$timetoEventData$eStatus == 1
exptime <- boxplot(rrAnalysisTrain$timetoEventData$expectedTime[isevent]~rrAnalysisTrain$timetoEventData$class[isevent],plot=FALSE)
obstime <- boxplot(rrAnalysisTrain$timetoEventData$eTime[isevent]~rrAnalysisTrain$timetoEventData$class[isevent],plot=FALSE)
classnames <- attr(rrAnalysisTrain$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
| 1Q |
168 |
92 |
286 |
176 |
| Median |
285 |
180 |
433 |
181 |
| 3Q |
431 |
329 |
447 |
271 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Cal. Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

Cox Calibration
op <- par(no.readonly = TRUE)
calprob <- CoxRiskCalibration(ml,lung,"status","time")
( 849.0437 , 781.9918 , 1.468556 , 121 , 140.6965 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
By Risk Categories
obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
| Total |
121 |
100.4 |
144.6 |
96.2 |
1.26 |
1.044 |
1.50 |
0.0142 |
| low |
83 |
66.1 |
102.9 |
64.2 |
1.29 |
1.029 |
1.60 |
0.0244 |
| 90% |
38 |
26.9 |
52.2 |
32.3 |
1.18 |
0.833 |
1.62 |
0.2912 |
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))
plot(obsexp$Expected,obsexp$Observed,
xlim=c(minx,maxx),
ylim=c(minx,maxx),
main="Cal. Expected vs Observed",
ylab="Observed",
xlab="Expected",
col=rainbow(nrow(obsexp)),
log="xy")
errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

Time to Event Analysis
isevent <- rrAnalysisTrain$timetoEventData$eStatus == 1
exptime <- boxplot(rrAnalysisTrain$timetoEventData$expectedTime[isevent]~rrAnalysisTrain$timetoEventData$class[isevent],plot=FALSE)
obstime <- boxplot(rrAnalysisTrain$timetoEventData$eTime[isevent]~rrAnalysisTrain$timetoEventData$class[isevent],plot=FALSE)
classnames <- attr(rrAnalysisTrain$timetoEventData,"ClassNames")
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
| 1Q |
168 |
92 |
230 |
142 |
| Median |
285 |
180 |
348 |
145 |
| 3Q |
431 |
329 |
359 |
218 |
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
xlab="Mean Expected Time",
ylab="Mean Observed",
xlim=c(minv,maxv),
ylim=c(minv,maxv),
main="Cal. Estimated Time to Event",col=rainbow(length(classnames)))
lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

Cross-Validation
rcv <- randomCV(theData=lung,
theOutcome = Surv(time,status)~1,
fittingFunction=BSWiMS.model,
trainFraction = 0.95,
repetitions=200,
classSamplingType = "Pro"
)
.[++].[+++].[+++].[+++].[++].[++].[+].[+++].[+++].[+++]10 Tested: 74
Avg. Selected: 3.5 Min Tests: 1 Max Tests: 4 Mean Tests: 1.351351 . MAD:
0.4864521
.[+++].[+++].[+].[++].[+++].[+++].[+].[++].[+++].[+]20 Tested: 118
Avg. Selected: 3.35 Min Tests: 1 Max Tests: 5 Mean Tests: 1.694915 .
MAD: 0.4734697
.[+++].[+++].[++].[++].[+].[++].[++].[+++].[++++].[+++]30 Tested: 144
Avg. Selected: 3.4 Min Tests: 1 Max Tests: 5 Mean Tests: 2.083333 . MAD:
0.4749689
.[+++].[+++].[+++].[+++].[++].[++].[++].[+++].[+++].[+++]40 Tested:
154 Avg. Selected: 3.475 Min Tests: 1 Max Tests: 7 Mean Tests: 2.597403
. MAD: 0.4785025
.[+++].[++].[+++].[+].[+++].[++].[+].[+++].[+++].[++]50 Tested: 159
Avg. Selected: 3.44 Min Tests: 1 Max Tests: 8 Mean Tests: 3.144654 .
MAD: 0.4783724
.[++].[++-].[++++].[+].[++].[+++].[+++].[++].[+++].[+++]60 Tested:
165 Avg. Selected: 3.45 Min Tests: 1 Max Tests: 9 Mean Tests: 3.636364 .
MAD: 0.4777576
.[++].[++].[+++].[++].[+].[+++].[+++].[+++].[++++].[++++]70 Tested:
165 Avg. Selected: 3.471429 Min Tests: 1 Max Tests: 9 Mean Tests:
4.242424 . MAD: 0.4783279
.[++++].[++++].[+++].[+++].[++-].[++++].[+++].[+++].[+++].[++-]80
Tested: 165 Avg. Selected: 3.55 Min Tests: 1 Max Tests: 10 Mean Tests:
4.848485 . MAD: 0.4782539
.[+++].[++].[+++].[+++].[+++].[++].[+].[+].[+].[+++]90 Tested: 166
Avg. Selected: 3.511111 Min Tests: 1 Max Tests: 11 Mean Tests: 5.421687
. MAD: 0.4781686
.[++++].[++].[+++].[+++].[++++].[+++].[++++].[+++].[++].[+]100
Tested: 167 Avg. Selected: 3.55 Min Tests: 1 Max Tests: 12 Mean Tests:
5.988024 . MAD: 0.4766993
.[++].[+].[+++].[+++].[+++].[+++].[+++].[++].[++++].[++]110 Tested:
168 Avg. Selected: 3.554545 Min Tests: 1 Max Tests: 14 Mean Tests:
6.547619 . MAD: 0.476821
.[++].[++].[++].[++].[+++].[++].[+].[+++].[+].[++]120 Tested: 168
Avg. Selected: 3.508333 Min Tests: 1 Max Tests: 14 Mean Tests: 7.142857
. MAD: 0.4766714
.[+++].[+++].[+++].[+++].[+++].[++].[++].[++].[+].[+++]130 Tested:
168 Avg. Selected: 3.507692 Min Tests: 1 Max Tests: 15 Mean Tests:
7.738095 . MAD: 0.4760375
.[+++].[++-].[++].[++].[+].[++].[+++].[+++].[++].[++]140 Tested: 168
Avg. Selected: 3.485714 Min Tests: 1 Max Tests: 16 Mean Tests: 8.333333
. MAD: 0.4760829
.[+++].[+].[+++].[+++].[+].[+++].[+++].[+++].[++].[+++]150 Tested:
168 Avg. Selected: 3.486667 Min Tests: 2 Max Tests: 16 Mean Tests:
8.928571 . MAD: 0.4764672
.[+++].[++].[++].[+].[+].[++].[+].[+++].[+++].[++]160 Tested: 168
Avg. Selected: 3.45625 Min Tests: 2 Max Tests: 19 Mean Tests: 9.52381 .
MAD: 0.476116
.[+++].[+++].[++].[+++].[++].[++].[++].[+++].[++].[+++]170 Tested:
168 Avg. Selected: 3.458824 Min Tests: 2 Max Tests: 20 Mean Tests:
10.11905 . MAD: 0.476003
.[++-].[+++].[++].[++].[+++].[+].[+++].[+++].[+++].[++++]180 Tested:
168 Avg. Selected: 3.466667 Min Tests: 2 Max Tests: 20 Mean Tests:
10.71429 . MAD: 0.4755427
.[+++].[++].[+].[+++].[++].[++].[+++].[++].[+++].[++]190 Tested: 168
Avg. Selected: 3.457895 Min Tests: 3 Max Tests: 21 Mean Tests: 11.30952
. MAD: 0.4754696
.[++].[++-].[++-].[++].[+++].[+++].[+++].[+++].[+++].[++-]200 Tested:
168 Avg. Selected: 3.46 Min Tests: 3 Max Tests: 23 Mean Tests: 11.90476
. MAD: 0.4756583
stp <- rcv$survTestPredictions
stp <- stp[!is.na(stp[,4]),]
bbx <- boxplot(unlist(stp[,1])~rownames(stp),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(unlist(stp[,2])~rownames(stp),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(unlist(stp[,4])~rownames(stp),plot=FALSE)$stats[3,],h0)
rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names
names(times) <- bbx$names
rrAnalysisTest <- RRPlot(rdatacv,atRate=c(0.90),
timetoEvent=times,
title="Test: Lung Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






Calibrating the test results
rdatacv <- cbind(status,prob,times)
calprob <- CalibrationProbPoissonRisk(rdatacv)
( 858.7891 , 781.9918 , 1.485412 , 121 , 143.3832 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
timeinterval <- calprob$timeInterval;
rdata <- cbind(status,calprob$prob)
rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=times,
title="Calibrated Test: Lung",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)





